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Abstract. The linear stability of a shear layer in a highly 
diffusive stably stratified atmosphere has been investi- 
gated. This study completes and extends previous works 
by Dudis (1974) and Jones (1977). We show that (i) an 
asymptotic regime is reached in the limit of large thermal 
diffusivity, (ii) there exists two different types of unstable 
modes. The slowest growing modes correspond to predom- 
inantly horizontal motions and the fast growing ones to 
isotropic motions, (iii) the physical interpretation of the 
stability properties in this asymptotic regime is simple. 
Applications to the dynamics of stellar radiative zone are 
discussed. 
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1. Introduction 

The standard model describing stellar evolution neglects 
the role of macroscopic motions in the radiative layers of 
stars. During the past twenty years, however, determina- 
tions of surface chemical abundances and helioseismology 
data have shown that some mechanical mixing must be 
invoked to account for the observations. Current mixing 
theories involve either motions induced by the rotation or 
gravity waves generated at the boundary with the convec- 
tive zones. 

Hydrodynamical instabilities due to differential rota- 
tion play an essential role in this context. Turbulent mo- 
tions arise through the non-linear development of hydro- 
dynamical instabilities and whenever present, they ensure 
an efficient mixing of chemical and angular momentum. 
The study of instabilities is therefore necessary to deter- 
mine the occurrence of turbulent mixing in stellar radia- 
tive zones. Moreover, in the absence of a better under- 
standing of turbulent motions, the properties of the insta- 
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bility are often used to infer the properties of the subse- 
quent turbulent motions. 

In a rotating fluid of uniform density, a hydrodynam- 
ical instability occurs if the specific angular momentum 
decreases outwards the rotation axis or if the rotational 
velocity varies along the rotation axis. Instability can also 
be triggered by a velocity shear, a typical example being 
the Kelvin-Helmholtz instability which develops at the in- 
terface between two streams of different velocities (see for 
example the reference books by Chandrasekhar 1961 and 
Drazin 1981). In principle, these different types of insta- 
bility can also take place inside stellar radiative zones. 
But there, the stable stratification will strongly oppose 
to them. Indeed, the time scale characterizing the action 
of the restoring buoyancy force is generally much smaller 
than the dynamical time scale (the inverse of the 

Brunt- Vaisala, frequency is of the order of 10 3 s whereas 
is equal to 10 5 5 s for the sun and is of the order of 
10 4 s for a fast rotating star with an equatorial velocity 
v cq = 200km s~ 1 and a radius R = 3Rq). 

This does not mean that instabilities are systemati- 
cally suppressed by stable stratification. But, as we shall 
see, this puts strong constraints on the type of motions 
which can arise from these instabilities. 

A possibility to avoid the action of the buoyancy force 
is to consider displacements limited to the surfaces of con- 
stant entropy. Such purely two-dimensional motions are 
not affected by the buoyancy force and might be ampli- 
fied if a differential rotation exists along these surfaces. 
However, even if the early development of the instability 
may not be affected by the stabilizing buoyancy force, the 
non-linear stage should produce three-dimensional mo- 
tions which in turn will be affected by the stable strat- 
ification. 

Another possibility is to consider three-dimensional 
perturbations whose vertical length scale is small enough 
to be significantly affected by the thermal diffusion. Ther- 
mal diffusion has indeed a destabilizing effect in stably 
stratified atmospheres as was first recognized in a geo- 
physical context by Townsend (1958). When a fluid par- 
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eel is displaced vertically from its hydrostatic equilibrium 
position, thermal diffusion reduces the temperature de- 
parture between the fluid parcel and its environment. For 
incompressible motions, the density departure and so the 
amplitude of the restoring buoyancy force are reduced by 
the same factor. Thermal diffusion has therefore a desta- 
bilizing effect because it weakens the stabilizing effect of 
the stable stratification. 

In the present paper, we investigate in details this 
destabilizing effect in the case of shear instabilities. Our 
work follows previous theoretical studies which we sum- 
marize now. Let's first recall that, in the absence of diffu- 
sive and viscous effects, there exists a criterion, derived in 
the context of the linear stability theory and valid for any 
parallel flows U = U(z)e x embedded in some vertically 
stably stratified Boussinesq atmosphere (z is the vertical 
coordinate and e x a horizontal unit vector) . This criterion 
established by Miles (1961) states that the flow is stable 
to infinitesimal perturbations if the Richardson number 
defined as, 



is everywhere larger than 1/4. 

Such a criterion does not exist in the presence of ther- 
mal diffusion. Existing results consist either in general 
criteria based on phenomcnological arguments or appli- 
cations of the linear stability theory to a particular flow 
i.e. a particular velocity profile within a particular atmo- 
sphere. 

The general criterion derived by Zahn (1974) states 
that the shear layer is stable if, 

1 

RiPe > -, (2) 

where Pe — uI/k is a Peclet number associated with 
an eddy of size I and velocity u. As expected, pertur- 
bations that would be stable in a perfect fluid because 
Ri > 1/4 can now be unstable provided Pe < l/(4i?i), 
i.e. if the thermal diffusion time scale of the eddy is small 
enough compared to its turnover time scale l/u. More re- 
cently, Maeder (1995) obtained a similar expression by 
extending to the diffusive case the energy considerations 
which allowed Chandrasekhar's derivation of Miles' crite- 
rion (Chandrasekhar 1961, see also Miles 1986). 

In the context of the linear stability theory, ther- 
mal diffusion effects have been studied independently 
by Dudis (1974) and Jones (1977). Both authors have 
chosen the hyperbolic-tangent velocity profile, U(z) = 
A[/tanh(2/L), whose instability properties are well- 
known in the unstratificd case (Drazin 1981). However 
they made different choices for the temperature profile 
characterizing the stable stratification. 

Dudis considered a hyperbolic-tangent profile for the 
temperature T(z) = T + ATtanh(z/L). He found in par- 



ticular that, for small enough Peclet numbers, the stabil- 
ity depends only on the product of the Richardson and 
Peclet number RiPe and not on Ri and Pe separately, 
this result being consistent with the criterion proposed by 
Zahn (1974) (here Ri = {NL) 2 /AU 2 and Pe = AUL/k). 
However, it should be noticed that the thermal stratifi- 
cation chosen by Dudis is not adapted to investigate the 
small Peclet number regime. Actually, in this regime, the 
diffusion of the hyperbolic-tangent temperature profile is 
much faster than the growth of the instability, and this 
invalidates a starting assumption of this type of insta- 
bility study, namely that the basic state is a stationary 
solution of the governing equations or, at least, can be 
considered as such during the development of the insta- 
bility. The minimum growth time of the Kelvin- Hclmholtz 
instability (i.e. the inverse of the maximum growth rate) 
is known to be ~ 5L/AU whereas the thermal diffusion 
time scale of the temperature profile is t K = L 2 /n. Thus, if 
the Peclet number Pe = AUL/k is smaller than unity, the 
basic temperature profile would have been diffused much 
before the instability had time to grow. 

To avoid this problem, one must consider a station- 
ary solution of the heat equation and, for a Boussinesq 
fluid, this implies a temperature profile varying linearly 
with altitude. Jones (1977) chose such a profile and then 
studied the stability in the range 1/4 < Ri < 1 as well 
as in the limit of very large horizontal wavelengths of the 
disturbance, k x , and very large Richardson number. The 
stability criterion in this double limit is k x PeRi > 0.086. 
However, the partial coverage of the parameter space (k x , 
Pe, Ri) did not allow to investigate the small Peclet num- 
ber regime, and in particular, the existence of the asymp- 
totic state suggested by Dudis (1974). 

Here we present a normal mode stability analysis of a 
hyperbolic-tangent shear layer embedded in a linearly sta- 
bly stratified Boussinesq atmosphere, with the objective 
of completing the work of Dudis (1974) and Jones (1977). 
First, we determine the conditions of marginal stability 
for a much larger range of Richardson numbers, Peclet 
numbers and Reynolds numbers. Then, we concentrate on 
the small Peclet number regime and investigate for the 
first time the growth rates of the unstable modes as well 
as their spatial properties. We also provide simple physi- 
cal interpretation of the stability properties by comparing 
the relative effects of shear, stable stratification and ther- 
mal diffusion. Lignieres (1999) recently derived an asymp- 
totic form of the Boussinesq equations in the limit of small 
Peclet number and showed that, in this asymptotic regime, 
thermal diffusion and stable stratification combine in a 
single process. We shall see that this property greatly sim- 
plifies the interpretation of our results. 

The paper is organized as follows. The equations gov- 
erning the evolution of perturbations are derived in Sect. 
3.1 and the numerical method to solve the corresponding 
eigenvalue problem is presented in Sect. 3.2. Results are 
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presented and interpreted in Sect. 4. In Sect. 5, they are 
summarized and discussed in an astrophysical context. 

2. The mathematical model 

2.1. Governing equations 

From now on, dimensional quantities are denoted with an 
overbar. We consider a parallel flow U = AU tanh(z)e x 
in an unbounded vertically stratified atmosphere. The z 
axis refers to the direction of gravity, while the x axis cor- 
responds to the stream direction. The equation of state is 
that of a Boussinesq fluid, g — q~q(1 — {3(f — To)) and the 
thermal stratification is given by the conductive profile 
T cq (z) — Tq + ATz/L. Here (3 is the thermal expansion 
coefficient and, go and To are references density and tem- 
perature. 

Using, L, AU and AT as units of length, velocity and 
temperature, respectively, the non-dimensional Boussi- 
nesq equations (Spiegel & Veronis 1960) read: 



1 9 

u • Vu = -Vp + Ride z H V 2 u 

Re 



du 

~dt 

30 1 

— + u • ve + w = — v 2 #, 

at Pe 
V-u = 0, 



(3) 

(4) 
(5) 



ve y + we z is the velocity vector, p the 



where, u = ue 
pressure and 8(x, y, z) = T(x, y, z) — T cq the temperature 
deviation from the temperature profile T eq (z). In the heat 
equation, the third term of the left hand side corresponds 
to the vertical advection of temperature against the mean 
temperature gradient dT cq {z) / dz which is equal to unity 
in our dimensionless units. 

The system is governed by three non-dimensional num- 
bers, the Richardson number, Ri, the Reynolds number, 
Re, and the Peclet number, Pe, respectively defined as 



Ri 



NL 
W 



Re 



AUL 



where, N = yj3gAT j X) is the Brunt- Vaisala frequency 
associated with the conductive temperature profile, v the 
molecular viscosity and R the thermal diffusivity. These 
non-dimensional numbers compare the dimensional time 
scales associated with the dynamics ts — L/AU, the 
stable stratification i N = 1/JV, the viscous dissipation 
t v = L 2 /D, and the thermal diffusion t K = L 2 /k. 

Due to the horizontal homogeneity, perturbations are 
resolved into modes of the form, 



Pe 



AUL 



U! 



{x,y,z) = ^(w{z)e i{ - k ^ +k y y - k ^), 



(6) 



where 3f() denotes the real part of a complex number, 
w(z) is the vertical profile of the mode (here the vertical 
velocity), k x and k y are the horizontal wavenumbers, and 
c = c r + ici is a complex number which real part c r is 



the phase velocity and which imaginary part Cj controls 
the growth or the decay of the mode. If a is positive, 7 = 
k x Ci is the growth rate of the unstable mode. Since two- 
dimensional perturbations contained in the plane (x, z) 
are more unstable than three-dimensional ones (Gage & 
Reid 1968), the spanwise velocity and the wavenumber 
k y can be taken equal to zero. Previous investigations on 
the stability of parallel shear flow indicate that the most 
unstable modes have symmetry properties according to 
the symmetry of the basic velocity profile. Here, due to the 
antisymmetry of the basic profiles with respect to z = 0, 
we can assume that w(z) = w*(—z) and <j>(z) = (f>*(—z), 
where cj>* is the complex conjugate of <fr. This implies that 
c r = 0. After linearization and some simple manipulations, 
one gets the two equations, 



ik x (U - c)M + ik x U"]w = Rikl$, 



[— - ik x (U- c)]$ = w, 
Pe 



(7) 



(8) 



where 4> is the vertical profile of temperature disturbance 
and M denotes the operator 4-^ — k x . The functions U and 
U" are the basic velocity profile and its second derivative, 
respectively. Together with the boundary conditions 



w,<p— >0 as z — > ±00, 



(9) 



Eqs. (§ and (§) fo rm an eigenvalue problem, denoted 
(El), which has been solved numerically. 

The eigenvalue problem depends on five param- 
eters k x , Ci, Ri, Pe, Re and may be written as 
T{k x ,Ci,Ri,Re,Pe) = 0. The hypersurface of marginal 
stability separating stable and unstable domains in the 
parameter space is given by the condition a = 0. 

2.2. Numerical model 

The linear system of ordinary differential equations ^ 
and (||) is solved with a numerical code that makes use 
of a finite-differences technique with deferred correction 
coupled to a Newton iteration. The step size can be dy- 
namically refined in the regions where the eigenfunctions 
present strong spatial gradients. 

Boundary conditions are applied at some distance z max 
from the shear layer. Then, we increase z m ax until the 
eigenvalues do no longer depend on the boundary condi- 
tions. Depending on the particular eigenfunction, z max w 
20 or 30 appeared to be sufficient to obtain an accurate 
solution. In the inviscid limit, the differential equations 
are singular at z — for neutral modes (cj = 0). Follow- 
ing Dudis (1974), we overcame this difficulty by solving 
the non-singular problem for Cj 7^ and then by decreas- 
ing Ci until the first two digits of the eigenvalue are not 
affected by the decay of Cj. These procedures have been 
successfully tested by recovering Dudis results. 
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3. Results 

We first present an overall picture of the stability prop- 
erties for a wide range of Richardson and Peclet num- 
bers, restricting ourselves to the inviscid case (subsection 
3.1). Then we concentrate on the regime of small Peclet 
number (subsection 3.2). The growth rates and the spatial 
structure of the unstable modes are studied in this regime 
(subsection 3.3) as well as the effect of a large Reynolds 
number (subsection 3.4). 

3.1. Overall presentation of the stability domains in the 
inviscid case 

The marginal stability curves are shown in Fig. 1 for 
a wide range of Peclet numbers: Pe = +oo, Pe = 1, 
Pe = 10 _1 , Pe = 10~ 2 and Pe = 10~ 3 . For any fixed value 
of Pe, the right hand side of the marginal curve corre- 
sponds to stable perturbations whereas the left hand side 
corresponds to unstable ones. The Pe = +oo curve is the 
analytical solution Ri c = k 2 (l — k 2 ) determined by Drazin 
(1958) (see also Chandrasekhar, 1961) in the non-diffusive 
limit. Each of the other curves has been drawn from about 
thirty numerically computed triplets (k x , Ri, Pe) between 
k x = 0.01 and k x = 1. 

It can be seen on Fig. 1 that a large thermal diffu- 
sivity has a destabilizing effect since the unstable domain 
increases as the Peclet number decreases. This is easily ex- 
plained as follows: When buoyancy is negligible, the tem- 
perature behaves like a passive scalar and thermal diffu- 
sion has no effect on the dynamics. Then, a first condition 
for a significative effect of thermal diffusion is that the flow 
be affected by the stable stratification, that is iff < Ig or 
equivalently Ri > 1. The second condition is that thermal 
diffusion reduces the stabilizing effect of the stratification. 
This will occur if thermal diffusion acts faster than sta- 
ble stratification, that is if I K < iff. Consequently the 
necessary conditions for a significant effect of the thermal 
diffusion are, 



t K < tjv < ts- 



(10) 



If one considers a shear layer characterized by a Richard- 
son number close to unity i.e. such that In ~ Ig, condition 
( |i"(i| ) states that thermal diffusion will reduce the stabiliz- 
ing effect of stable stratification if I K < Iff — ig that is if 
Pe < 1. Now, if the shear layer is very strongly stratified 
that is Ri >■ 1, or Iff <C Is, condition ( |Io| ) requires that 
t K <C is or Pe <C 1 to destabilize the layer. These sim- 
ple arguments account for the mean displacement of the 
neutral curves towards larger Richardson number, start- 
ing from the Pe — 1 curve. Thus, condition ([l0|) can be 
said to explain the global behaviour of the stability picture 
shown in Fig. 1. 

However, another point apparent on Fig. 1 is that large 
horizontal wavelengths are stable for Pe = +oo but be- 
come unstable as Pe decreases. They become even the 



most unstable modes at fixed values of Pel This property 
is harder to understand because thermal diffusion acts pri- 
marily on small scale perturbations. 

We give here a partial answer by showing how pertur- 
bations with large horizontal wavelength can be strongly 
affected by thermal diffusion: Again, in order to weaken 
the stable stratification, thermal exchanges must act be- 
fore stable stratification. This condition has been previ- 
ously expressed by i K < In but this inequality does not 
take into account the anisotropy of the buoyancy force. 
Actually, the effect of stable stratification tends to van- 
ish when motions become more and more horizontal. By 
contrast, thermal diffusion remains efficient for these mo- 
tions since predominantly horizontal motions must vary 
on a relatively small vertical scale to ensure mass conser- 
vation (V-u = 0). Therefore, thermal diffusion tends to be 
more efficient than stable stratification for predominantly 
horizontal motions. 

This reasoning can be formalized by considering in- 
finitesimal perturbations of the form w oc exp[i(k x x+k z z)] 
in a linearly stratified atmosphere at rest. The stable strat- 
ification time scale then corresponds to the inverse of grav- 
ity waves frequency, 



t N {k x , k z ) = yjkl + k 2 Jk x VRi, (11) 

whereas the thermal diffusion time scale is simply, 

t K (k x ,k z ) = Pe/(k 2 x + k 2 z ). (12) 

Comparing both time scales shows that, whatever the 
value of the Peclet number, thermal diffusion acts much 
faster than stable stratification when | k x /k z \ — \ w/u \ 
goes to zero, that is when motions are predominantly hor- 
izontal. 

To show that the motions associated with the marginal 
modes at large horizontal wavelengths are in this case, we 
calculated the square root of the ratio between the vertical 
and horizontal kinetic energy of the modes. The ratio is 
given by E w /E u where, 



+ 00 r2-K/k x 



E„ 



r+OO ri 
J-oo JO 



dxdz 



r-+oo rSir/kx, 9 2\ 
J-oo Jo ( W + U ) 



dxdz 



(13) 



and E u is defined accordingly. It has been determined for 
marginal modes taken along the Pe = 1 neutral curve and 
the results are presented on Fig. 2. It clearly appears that 
disturbances with large horizontal wavelengths are associ- 
ated with strongly anisotropic, predominantly horizontal 
motions. Actually, in the limit of small k x , the eigenfunc- 
tions w and <f> do not vary anymore so that the vertical 
scale height of the disturbance remains fixed while its hor- 
izontal wavelength increases. This explains why thermal 
diffusion can affect their stability, regardless of the value 
of the Peclet number. What is still unclear is why they 
are the most unstable disturbances. We shall come back 
to this point at the end of the next subsection. 




3.2. The small Peclet number regime 

In the following, we concentrate our analysis on the small 
Peclet number regime which is most relevant for stel- 
lar astrophysics. As suggested by Dudis' work, we plot 
once again the neutral curves corresponding to Pe = 
1, 10 _1 , 10~ 2 , 10~ 3 but this time as a function of RiPe. 
The outcome is displayed on Fig. 3. We observe that three 
of the marginal stability curves nearly collapse into a sin- 
gle curve. While differences between the Pe = 10 -2 and 
Pe = 10~ 3 curves are not distinguishable, the Pe = 1CP 1 
curve slightly deviates from those. 

This remarkable property strongly suggests the exis- 
tence of an asymptotic state for small Peclet numbers 
in which the governing equations only depend on the 
non-dimensional number, RiPe. Dudis (1974) had already 
found that the eigenvalue problem for neutral modes de- 
pends only on the product RiPe if the mean flow advec- 
tion is neglected in Eq. (j^). Using asymptotic expansions 
and starting from the Boussinesq equations, Lignieres 
(1999) recently obtained a general form of the asymptotic 
equations in the limit of small Peclet numbers. It is shown 



that, if u and 6 behave like Taylor series for small Peclet 
numbers (u = Uo + Peu± + ...,8 = 9q + Pe6\ + ...), the 
solutions of the Boussinesq equations are identical to the 
solutions of the following system up to the first order in 
Pe: 



du 



1 



— + u • Vu = -Vp + Ripe z + — V z u, 



Of 

V 2 ip = w 
V ■ u = 



Re 



(14) 

(15) 
(16) 



where rp = 9/ Pe are the temperature deviations scaled by 
Pe and 



R = RiPe = 



tst K 



(17) 



In this context, the linear equations governing the shear 
layer stability become: 

M 2 

[ ikJU - c)M + ik x U"}w = Rklip, (18) 

Re 
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Fig. 3. The marginal stability curves shown on Fig. 1 are now plotted as a function of RiPe. The collapse of the 
Pe = 10~ 2 and Pe = 10~ 3 curves suggests the existence of an asymptotic regime in the limit of small Peclet number. 
Moreover, all the curves tends towards the curve k x RiPe ~ 0.117 as k x goes to zero, suggesting the existence of 
another asymptotic regime in the limit of small horizontal wavenumber 



M^ = w, (19) 

which, as expected, only depend on the non-dimensional 
numbers, R = RiPe and Re. The eigenvalue problem 
formed by the above equations and the boundary con- 
ditions fl9j) has been solved numerically for neutral (c = 0) 
and inviscid (Re = +00) modes. We found that the re- 
sulting marginal stability curve can not be distinguished 
from the Pe = 1CP 2 and Pe = 1CP 3 curves obtained in 
the context of the Boussinesq equations. 

The existence of this asymptotic regime has various in- 
terests for the stability analysis. First, we can concentrate 
on the asymptotic limit and do not have to repeat the 
analysis for different values of the Peclet number. Second, 
it turns out that the stability is no longer governed by 



three processes (shear, stable stratification and thermal 
diffusion) but only by two, the shear destabilization and a 
stablizing process combining the stable stratification and 
thermal diffusion effects. 

To understand why, let us consider the differences be- 
tween the heat conservation of the Boussinesq equations 
(Eq. (Q)) and the simplified heat balance of the small- 
Peclet- number approximation (Eq. (JTsj) ) . This simplified 
balance occurs because, unlike other advection terms, the 
vertical advection against the mean stratification does not 
vanish when the high thermal diffusion strongly reduces 
the temperature deviations. Indeed, as fluid parcels go up 
(or down) in a mean temperature gradient, the ampli- 
tude of the associated temperature deviations tends to in- 
crease continuously. While this process acts as a source of 
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Fig. 2. Square root of the ratio between the vertical and 
horizontal kinetic energy of neutral modes taken along the 
Pe = 1 marginal stability curve. This figure shows that 
the motions associated with these modes become more 
and more horizontal as k x goes to zero 



temperature deviations, thermal diffusion tends to smooth 
them out. Eq. (|l5| ) represents the balance between both 
processes. In this state, vertical advection against mean 
temperature gradient and thermal diffusion no longer act 
individually to determine the amplitude of the buoyancy 
force; the process which combines their action will be 
called the small-Peclet-number buoyancy in the following 
of the paper. Its elementary properties have been ana- 
lyzed in Lignieres (1999) and are briefly summarized here. 
It is purely dissipative in the sense that the work of the 
buoyancy force integrated over the whole domain is always 
negative, provided temperature deviations vanish on the 
outer boundaries. Its time scale is Ib = tj^/i K , although 
the process being strongly anisotropic, the time scale of 
the small-Peclet-number buoyancy is better described by, 



tB{k x ,k z ) = 



(hi 



Rkl 



(20) 



which corresponds to the time scale necessary to dissipate 
an infinitesimal perturbation of the form w oc exp[i(k x x + 
k z z)\ in the atmosphere at rest. 

Coming back to the stability analysis, we shall now 
compare the stabilizing effect of the small-Peclet-number 
buoyancy with the destabilizing one of the shear. 

Let us consider the instability of large horizontal wave- 
length disturbances. In the previous subsection, we have 
shown that the vertical length scale of these disturbances 
remains fixed as their horizontal wavenumber goes to 
zero. Then, according to expression (p0|), the small-Peclet- 
number buoyancy time scale increases like l/k x in the 
same limit. By comparison, the time scale of the insta- 
bility (the inverse of the growth rate in a fluid of constant 
density) is known to be proportional to l/k x for small 



k x (Drazin 1981). Consequently, the stabilizing effect de- 
creases faster than the destabilizing one as k x vanishes. 
This explains why predominantly horizontal perturbations 
are unstable and why they are increasingly unstable as k x 
tends to zero. 

It is worth mentioning the difference with the non dif- 
fusive case (Pe = +oo). In this case also, the stable strat- 
ification is inefficient for predominantly horizontal distur- 
bances. But, by contrast with the diffusive situation, the 
stable stratification time scale increases as l/k x (see Eq. 
(|TT|)), that is as fast as the instability time scale. The re- 
sult is that disturbances with large horizontal wavelength 
are stable (sec the Pe = +oo neutral-curve on Fig. 1). We 
therefore conclude that the shapes of the non-diffusive and 
highly diffusive marginal stability curves differ because the 
non-diffusive buoyancy has a stronger effect on predom- 
inantly horizontal motions than the small-Peclet-number 
buoyancy. In other words, the action of the small-Peclet- 
number buoyancy is more anisotropic than that of the 
non-diffusive buoyancy, and this triggers the instability of 
predominantly horizontal motions. 

The small-Peclet-number buoyancy allows also to un- 
derstand why, as indicated in Fig. 3, the marginal stabil- 
ity curves converge towards k x R = 0.117 when k x goes 
to zero. Since stabilizing and destabilizing effects balance 
each other for marginal modes, the marginal curve can be 
characterized by an equality between the dynamical time 
scale and the small-Peclet-number buoyancy time scale. 
In the limit of vanishing k x , this reads: 



1 

k x 
that is, 



BR 1 



k x R = constant, 



(21) 



(22) 



since the vertical scale height of these modes remains fixed. 

Although outside the regime of small Peclet numbers, 
we observe that the Pe = 1 marginal curve also converges 
towards k x R = 0.117 for small k x . This is due to the fact 
that modes associated with predominantly horizontal mo- 
tions satisfy the simplified heat equation ( [i"5| ) whatever 
the value of the Peclet number. In the limit of small k x , 
the advective term UdO/dx vanishes naturally while the 
temporal variation 88/ 'dt is also negligible because, as we 
have seen before, thermal diffusion dominates stable strat- 
ification for any value of the Peclet number provided the 
motions are sufficiently horizontal. 

From a mathematical point of view, the convergence 
towards k x R — 0.117 means that the eigenvalue problem 
has reached an asymptotic regime in the limit of small k x . 
Effectively, when terms proportional to k x are neglected 
in Eqs. (18) and (p9|), the asymptotic eigenvalue problem 
depends only on the parameter k x R. Jones (1977) had al- 
ready obtained this asymptotic eigenvalue problem, start- 
ing from Eqs. ((?]) and (||) and considering the double limit 
k x — > and Ri — > +oo. He found a slightly smaller value, 
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Fig. 4. Growth rate isocontours in the limit of small Peclet 
numbers. The bold curve corresponds to marginal stabil- 
ity (7 = 0). Note that the maximum growth rate in the 
unstratified case is 70 — 0.189 

k x R — 0.086, most probably because of the difference in 
the boundary conditions. In his study, the requirement 
that perturbations vanish at z = ±4 decreases a little bit 
the domain of instability. 

3.3. Growth rates in the small Peclet number regime 

In the unstratified case, it is well known that the growth 
rate of the instability is maximum for a particular mode 
(70 ~ 0.189 at k x ~ 0.447, see for example Drazin 1981). 
In the same way, one would like to know whether, for a 
fixed value of R = RiPe, a particular mode will grow 
faster than the others. In that relevant issue will 

be to quantify the reduction of the maximum growth rate 
due to the small-Peclet-number buoyancy. 

Contours of constant growth rate have been deter- 
mined inside the unstable domain and are presented on 
Fig. 4. Unsurprisingly, the growth rates decrease when 
RiPe increases, that is when the stabilizing effect of the 
buoyancy force becomes stronger. For example at RiPe ~ 
5, the maximum growth rate is ~ 0.001, that is smaller 
by more than two order of magnitude than the maximum 
growth rate of the unstratified case. 

It can be inferred by looking at Fig. 4 that a maxi- 
mum growth rate exists for each fixed value of R = RiPe. 
This maximum value denoted 7 max is shown on Fig. 5 as 
a function of R. It appears that a rapid decrease of the 
maximum growth rate is followed by a much slower one 
after an abrupt change of slope at R ~ 0.253. This change 
of behaviour is due to the existence of two different types 
of modes as supported by Fig. 6 where, &™ ax , the hori- 
zontal wave number of the most amplified disturbance is 
plotted against 7 max . 
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Fig. 5. Maximum growth rate as a function of R = RiPe. 
Only the range < R < 1 is represented here 
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Fig. 6. Horizontal wavenumber of the most unstable 
modes 



Figs. 5 and 6 show that a first type of mode character- 
ized by horizontal wavenumbers close to 0.5 is dominant 
in the regime R < 0.253 whereas a second type of mode 
with much smaller wavenumbers becomes dominant once 
R > 0.253. The "jump" between both types of modes oc- 
curs at R ~ 0.253 when the curve 7 = f(k x ) has two equal 
local maxima. 

The vertical structure of both types of modes has been 
analyzed in order to characterize their geometrical prop- 
erties. For the modes which dominate at small values of 
R, the vertical decay outside the shear layer is quite rapid 
since their vertical scale height varies between L and 2L. 
These modes can be said to be isotropic as their vertical 
and horizontal length scale arc of the same order. 

By contrast, the modes which dominate for large val- 
ues of R are strongly anisotropic and become more and 
more so as R increases. Effectively, we found that their 
vertical scale height remains fixed to ~ 5.56L while their 
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00.0 



Fig. 7. Maximum growth rate as a function of R — RiPe 
in log-log coordinates 

horizontal wavelength ranges from 24L to +oo as R in- 
creases. Another property of this type of modes is that 
their growth rate and horizontal wavenumbers are related 
to R through scaling laws. This can be seen on Fig. 7 
where the relation between 7 max and R is plotted again 
using log-log coordinates and including large values of R. 
The scaling law, 



7n 



5. x lO'^R- 



(23) 



is rapidly established as R increases above 0.253. In the 
same way, it is found that the wavenumber of the most 
amplified mode and R are related by 



4.86 x 1Q- Z R- 



(24) 



These scaling laws are due to the fact that the eigenvalue 
problem reaches an asymptotic regime for small horizontal 
wavenumbers. Indeed, in the limit k x — > 0, the eigenvalue 
problem for non-neutral modes (k x Ci = 7 7^ 0) is gov- 
erned by two parameters only, a = k x R and b = j/k x . 
Relations (^3|) and fl24| ) then show that the most unstable 
mode corresponds to a particular solution of the asymp- 
totic eigenvalue problem. 

3-4- Effect of a large Reynolds number in the small 
Peclet number regime 

For a constant density fluid, the Kelvin- Helmholtz insta- 
bility is said to be inviscid because viscosity does not play 
any role in the mechanism which triggers the instability. 
Viscous dissipation can still influence the stability but, at 
least for the hyperbolic-tangent shear layer, its effect is 
rapidly negligible as the Reynolds number increases. This 
is not the case in a highly diffusive atmosphere since we 
found that viscous dissipation has a significant influence 
even if the Reynolds number is very large. This point is 
illustrated by Fig. 8, where the inviscid neutral curve is 
compared with the neutral curve calculated at Re — 1000. 



STABLE 




Fig. 8. Curves of marginal stability for Re = 1000 and 
Re = +00 in the small-Peclet-number regime 



The figure shows that viscosity essentially comes into play 
by stabilizing the large scale predominantly horizontal dis- 
turbances, its influence on the other type of perturbations 
being much weaker. 

This can be explained as follows: for wavelengths of 
the order of the shear layer thickness, the viscous time 
scale is much larger than the other time scales. However, 
as k x goes to zero, both the time scale of the small-Peclet- 
number buoyancy and the time scale of the instability be- 
come infinite while the viscous time scale tends towards a 
constant value, 



tj/ H ™ -^"^G 



(25) 



where H z denotes the vertical length scale of the predom- 
inantly horizontal modes. Then, at some point, viscous 
dissipation becomes dominant and stabilizes the pertur- 
bations. 

The critical values of R and k x have been determined 
for two other Reynolds numbers, Re = 1100 and Re = 
1500. They staisfy the following relations, 



Rait 

Re 



1.67 x 10~ 3 Rett 



crit 



19, 



so that the stability criterion is, 



RiPr 



R 
Re~ 



> 1.67 x 10" 



(26) 



(27) 



for the three Reynolds numbers considered. 

There are two different ways of recovering the rela- 
tions ( p6[ ) by simple arguments. First, by saying that all 
perturbations are stable when viscous damping balances 
the growth of the most unstable mode in the inviscid case. 
Equalizing the viscous time scale t v with the inverse of the 
maximum growth rate l/7 ma x yields R/Re — constant. 
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The second scaling law (g3[) shows that the corresponding 
wavenumber satisfies k x Re — constant. 

The same relations can be obtained without knowing 
the expression of the maximum growth rate. We first de- 
fine a stabilizing time scale which includes the effects of 
the small-Peclet-number buoyancy and the viscous dissi- 
pation as, 

(^k x ) 



It corresponds to the damping time of an infinitesimal 
perturbation proportional to exp[i(k x x + k z z)] in a qui- 
escent atmosphere described by the small-Peclet-number 
approximation. Then, the marginal stability curve can 
be characterized by an equality between this time scale 
and the instability time scale. In the limit of small k x , it 
turns out that the function R{k x ) possesses a maximum 
-Rmax = -Refc^/4, reached at k x Re — k 2 z . The above rela- 
tions follow since the vertical scale height of these modes 
remains fixed. 

For smaller Reynolds numbers (25 < Re < 150), Jones 
(1977) found a criterion similar to (p7|), although the con- 
stant is somewhat higher RiPr > 7x 10 -3 . This difference 
can be attributed to the boundary conditions. 

4. Summary and discussion 

The linear stability of a hyperbolic-tangent shear layer in 
a stably stratified atmosphere has been investigated for 
a wide range of Richardson, Peclet and Reynolds num- 
bers. Emphasis has been put on the regime of small Peclet 
numbers (i.e. high thermal diffusivity) relevant for stellar 
interior dynamics. Note that this regime is also impor- 
tant for the upper atmosphere of the Earth and Venus as 
Townsend (1958) and Dudis (1974) pointed out. 

In the inviscid case, we found that the shear layer is 
always unstable provided the horizontal wavelength of the 
perturbations is not bounded. It should be specified that 
the growth rate goes to zero as the horizontal wavelength 
goes to infinity. Viscous dissipation stabilizes this type of 
perturbation and so introduces a critical value of R = 
RiPe above which the shear layer is stable. 

The analysis of the growth rates in the inviscid small 
Peclet number regime revealed the existence of two differ- 
ent types of modes. A first type, isotropic and relatively 
fast growing, is dominant in the range < R — RiPe < 
0.253. By relatively fast, we mean that the growth rate is 
always larger than O.I70, where 70 is the maximum growth 
rate in the unstratified case. For larger values of R, the 
most unstable modes are predominantly horizontal. Their 
growth rates continue to decrease, following the scaling 
law 7 max = 5. x 1Q~ S Ri~ t~*, in dimensional units. 

Perhaps, the most striking result is the existence of 
an asymptotic regime where the mathematical and physi- 
cal analysis of the stability is considerably simplified. The 



physical simplification comes from the fact that the buoy- 
ancy force is no longer determined by two independent 
processes, temperature advection in the stratified atmo- 
sphere and thermal diffusion, but rather by a unique pro- 
cess combining both. The effect of this so-called small- 
Peclet-number buoyancy turns out to be purely dissipative 
and strongly anisotropic. Comparing its stabilizing effect 
to the destabilizing Kevin- Helmholtz mechanism leads to 
a straightforward interpretation of the stability properties. 

It should be noted that the small-Peclet-number buoy- 
ancy characterizes the response of a highly diffusive stably 
stratified atmosphere to any sufficiently slow motions. As 
such, it can potentially be applied to other types of mo- 
tions than those driven by a parallel shear flow. This in- 
cludes for example other types of double diffusive instabil- 
ities like the Goldreich-Schubert instability or the thermo- 
halinc instability. We note in particular that the most un- 
stable modes of the latter correspond to long thin columns. 
This is fully compatible with the properties of the small- 
Peclet-number buoyancy since the associated stabilizing 
time scale (^0|) becomes infinite for this type of motions 
(k x — > +00 and k z remaining fixed). Work is in progress 
to assess the usefulness of the small-Peclet-approximation 
in this context. 

Strictly speaking, the results presented here are only 
valid for the particular flow considered and assuming that 
the perturbations are infinitesimal. We expect neverthe- 
less that similar results would be obtained for any shear 
layer subject to the Kelvin-Helmholtz instability, that is 
when the vorticity profile possesses an extremum (Drazin, 
1981). But, in any case, our results can not be directly 
extended to the other velocity profiles and more impor- 
tantly if the basic flow is perturbed by finite amplitude 
disturbances. Effectively, since Reynolds' experiments on 
flows in a pipe (Reynolds 1883), it is well known that finite 
amplitude perturbations are able to trigger shear layer in- 
stabilities when infinitesimal ones can not. In this context 
the criterion we obtained appears as a sufficient - but not 
necessary - condition for instability. 

In principle, it can be applied as such, to shear layers 
inside stellar radiative zones. However, taking typical so- 
lar values for the thermal diffusivity and the Brunt- Vaisala, 
frequency (k = 10 7 cm 2 s _1 and N = 10~ 3 s _1 ), it is found 
that shear layers unstable according to this criterion would 
have vertical extent much smaller than the vertical reso- 
lution of helioseismology data (L/Rq < 10~ 5 ). To give 
an example, if one assumes that the shear across the solar 
tachocline can be modeled by a hyperbolic-tangent pro- 
file entirely embedded in the radiative zone, the associ- 
ated Richardson number would be about 2.5 x 10 3 and 
the Peclet number about 10 6 (we took 5 percent of the so- 
lar radius for the tachocline thickness and Ail = 10 _6 s _1 
for the differential rotation across the tachocline). Then, 
our criterion would predict stability because the small- 
est unstable horizontal wavelength according to k x = 
0.117 / RiPe is much larger than the solar radius. Note 
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again that shears across the tachocline could be subjected 
to finite amplitude instabilities as discussed in Michaud & 
Zahn (1998). 

For the time being, one relies on the criterion estab- 
lished by Zahn (1974) to decide about the non-linear sta- 
bility of a shear layer. This criterion is not in contradiction 
with our results. But this consistency check is not suffi- 
cient to prove its validity. Interestingly enough, the ener- 
getic consideration made by Maeder (1995) have shown 
that the criterion is consistent energetically. However, in 
the derivation by Zahn (1974) and by Maeder (1995), 
perturbations with a velocity scale u and a length scale 
I are introduced ab initio and it is not known whether 
such perturbations correspond to possible solutions of the 
equations of motion. Therefore, this criterion can be in- 
terpreted as a necessary condition for instability. 

Further investigations are clearly needed to better con- 
strain this criterion. In a geophysical context, numerical 
simulations of stably stratified sheared turbulence have 
been able to find critical Richardson numbers for which 
turbulence neither grows nor decays (see a review by 
Schumann, 1996). Such simulations should be performed 
at smaller Peclet numbers to approach stellar conditions 
and, as pointed out by Lignieres (1999), the small-Peclet- 
number approximation would be very useful in this con- 
text. Another possible approach is to assume that finite 
amplitude disturbances provoke a defect of the basic pro- 
file and then to analyze the linear stability of the per- 
turbed velocity profile as a function of the defect ampli- 
tude. This method has been used by Dubrulle & Zahn 
(1991) for the unbounded Couette flow in a constant den- 
sity fluid and could be extended to the stably stratified 
case at small Peclet number. 

As quoted in the introduction, beyond the stability of 
shear layer, the final objective is to estimate the angu- 
lar momentum and chemical element transport driven by 
shear layer turbulence. Assuming that small scale turbu- 
lence behaves like a diffusion process, Zahn (1992) pro- 
posed an expression for the turbulent viscosity on the ba- 
sis of arguments used to derive the stability criterion. This 
expression reads, 



u t = constant x k\ — — oc nRi a , 
N dr 1 



(29) 



where Ri g is the Richardson number formed with the lo- 
cal velocity gradient fdVl / df and the local Brunt- Vaisala 
frequency. Clearly our linear stability analysis gives no 
clues about the non-linear processes governing turbulent 
transport. Nevertheless, as the instability mechanism con- 
trols the energy injection into turbulence, it is natural 
to suppose that changes in the stability properties have 
a direct impact on turbulent transport. On dimensional 
grounds, we can form a vertical transport coefficient using 
the growth rate and the vertical length scale characteriz- 
ing the most unstable mode is 7 max -fff ■ For large values of 
R = RiPe (R > 0.253), 7 max follows the scaling law ( f23| ) 



whereas the typical vertical length of the most unstable 
mode remains fixed, so that Tmax-ff? oc RRi^ 1 . This ex- 
actly corresponds to the scaling of the turbulent viscosity 
(p9|). Such a close agreement is not found in the regime 
of small RiPe. However, the variation of the maximum 
growth with RiPe and the expression of the turbulent vis- 
cosity are still compatible. Actually, the main difference 
comes from the very existence of two different regimes 
which are not contained in the expression of the turbulent 
viscosity (29). One corresponds to rapid isotropic motions 
on the length scale of the shear layer, the other to slow 
nearly horizontal motions. 
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